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Abstract. The one-dimensional triplet contact process with diffusion (TCPD) model 
has been studied using fast multispin GPU Monte Carlo simulations. In particular, 
the particle density p and the density of pairs of neighboring particles p p have been 
monitored as a function of time. Mean field predictions for the time evolution of these 
observables in the critical point are p ~ t~ s and p p ~ t~ Sp with (5 = 1/3 and 5 P = 2/3. 
We observe that in the vicinity of the critical point of the model, the ratio p p /p tends 
to a constant, which shows that the one-dimensional TCPD model is not described by 
mean field behavior. Furthermore, our long simulations allow us to conclude that the 
mean field prediction of the exponent <5 is almost certainly not correct either. Since the 
crossover to the critical regime is extremely slow for the TCPD model, we are unable 
to pinpoint a precise value for 5, though we find as an upper bound 6 < 0.32. 
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1. Introduction 

The triplet contact process with diffusion (TCPD) model belongs to a set of closely 
related models of "fermionic" particles on a lattice that follow very simple dynamical 
rules. With fermionic we mean in this context that only one particle can be 
present at each site. These models have been studied extensively, because they were 
conjectured by Grassberger [TJ and Janssen [2] to belong to a relatively small number of 
dynamical universality classes, determined by coarse features such as the dimensionality, 
symmetries in the model, and conservation laws; in close analogy to universality classes 
in equilibrium statistical physics. 

Historically, the first of these models to be studied extensively is the Directed 
Percolation (DP) model. This model lends itself extremely well for computer 
simulations, and its exponents are therefore known with high accuracy. For example, 
using exact enumeration techniques, the exponent S was found to be 0.159464(6) [3], 
where 5 is defined through the time dependence of the particle density p of the system, 
starting from a state with a uniform high density: 

p~t~ s . (1) 

Even though this model is very easy from a numerical point of view, there are 
no theoretical predictions for these exponents, not even in one dimension. Since the 
Grassberger-Janssen conjecture says that the DP model belongs to a larger class of 
models with the same critical exponents, much effort has been undertaken to verify 
numerically the exponents in models that are also expected to be in the same universality 
class. 

In one interpretation of the DP model, particles are placed on a lattice, and follow 
two reactions: with a statistical rate p, each particle annihilates, and with a statistical 
rate 1 — p, each particle creates a new particle on an adjacent site, provided that it is 
vacant. Sometimes, the particles can also hop to neighboring lattice sites with a diffusion 
rate d. An extension of the DP model that has been studied thoroughly in literature 
is the pair contact process with diffusion (PCPD) model, in which the annihilation and 
procreation reactions can only take place if two particles are placed next to each other 
(and with d > 0). There has been a lot of discussion about the critical exponents, 
mainly 5, with estimates ranging from the DP- value 5 = 0.159 [H |5] to 8 — 0.35 [6]. 
In Ref. [7], we performed simulations of the PCPD model in which the triple product 
V = N x L x T, where N independent simulations are performed with a lattice of size L 
over a time range T, is much larger than in previous studies. The resulting data showed 
that in the case of the PCPD model, finite-time corrections are particularly severe, and 
with careful analysis we found more evidence suggesting a value close to the DP value, 
than the contrary. 

Here, we simulate the TCPD model, which makes annihilation and procreation 
conditional on triplets of particles. Currently, the claim in literature is that TCPD is 
fundamentally different from both the PCPD and DP models, as evidenced by different 
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values for exponents, e.g. 5 = 0.27(1) [8] and 5 = 1/3 [HI [TO] . We use the same efficient 
algorithm, described in [TTj . to reach a triple product (N x L x T) which is almost 
three orders of magnitude larger than in these previous studies. We find that for shorter 
times, the exponent 5 is very close to its mean field prediction of 1/3, but for very 
long times, this exponent starts to drift. We attribute this to finite-time effects that 
are even stronger than in the case of PCPD. Even with our fast algorithm and long 
simulations, we were not able to reach time scales that make it possible to retrieve an 
accurate estimate for 5. Importantly, however, we find that the DP value cannot be 
excluded, meaning that it is not yet proven that the TCPD model falls outside the DP 
universality class. 

2. Method 

Several slightly different versions of the TCPD model have been studied. We use a 
version with the reactions 3A — y AA and 3^4 — y 0, which is described by the following 
reactions and rates: 

AAAO -> AAAA 
OAAA -> AAAA 
AAA -»■ 000 
AO <-> OA 

In a straightforward implementation of this model, first the type of reaction is 
selected, based on a random number r e [0, 1), and depending on its value, one of the 
reactions is proposed (but not always carried out): 

• if r < d, a random pair of neighboring sites {i,i + 1} is selected; in case one site 
is occupied and the other is empty, the particle hops from the occupied site to the 
empty one. 

• else if r — d < p(l — d), a random triplet of sites {i, i + 1, i + 2} is selected; if all 
three sites are occupied, they are all made vacant. 

• else if r — d — p(l — d) < (1 — p)(l — d)/2, a random triplet of sites {i, i + 1, % + 2} is 
selected; if all three sites are occupied and site % + 3 is vacant, a particle is placed 
on this last site. 

• else, a random triplet of sites {i, i+2} is selected; if all three sites are occupied 
and site i — 1 is vacant, a particle is placed on this last site. 

Irrespective the type of reaction and its success, the time scale is incremented by 
At = 1/N. These steps are then iterated many times. As noted in the introduction, for 
the results presented here, we used as a basis an algorithm that leverages the power of 
graphics processing units (GPUs), which is described in detail in [TT] . 

Before one can actually measure the exponent 6, a prerequisite is to determine 
the critical annihilation rate p c . For low annihilation rates (p < p c ), the system will, 



each with rate 

with rate 
with rate 



-p)(l-d)/2 

p(l-d) 
d 



(2) 
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given enough time, settle with extremely high probability for a more or less constant 
particle density This regime is called the active regime. On the other hand, in the 
inactive regime with high annihilation rates (p > p c ), the particles will quickly die out. 
In between, exactly at p = p c , there is the critical regime, where the density decreases 
slower than in the inactive regime, following a power-law decay with a critical exponent 
5. 

Thus, it is necessary to first identify an estimate for the critical point p c , before we 
can estimate the critical exponent 5. To this effect we use as our main tool the effective 
exponent 8 e g function of time, defined as: 

log(p(ti)) - log(p(t 2 )) 



5eff{t = VUh) = : 7— : 7— , (3) 

log(ti) - log(t 2 ) 

and a similar expression for the effective exponent S e s iP for the pair density. 
Substituting the asymptotic behavior as given in eq. (pQ), we retrieve that 5 c q 5 
as t goes to infinity. The procedure is equivalent to numerical differentiation of 
<91og(p)/<91og(t). Equation (j3J) shows that there is still freedom in choosing t\ and 
t 2 . The trade-off is as follows: if we choose t\ closer to t 2 , the plot is generally more 
accurate, in the sense that features present in the analytical 5 e s curve are less smoothed 
out and lost that way. On the other hand, the curve is much more noisy. We found that 
in our case choosing t\jti ~ exp(3) gives good results. 

We find an estimate for p c by a manual binary search, which sounds more 
cumbersome than it is, because far from p c the 5 e s plots are very clearly recognizable 
as either sub- or super-critical. Closer to p c it gets much harder to estimate p c , because 
of the drift in the effective exponent, which gives our estimate for p c a larger error bar. 
Apart from the density p of the system, we can also measure the pair density p p . This is 
especially useful for testing the validity of mean field theory, as it predicts that p p = p 2 , 
and thus also S e s tP = 2S e s. 

For each value of p we performed at least iV = 200 independent simulations of runs 
up to T = 8 - 10 8 , with a lattice size of L = 2 21 = 2097152. We simulated at annihilation 
rates between p = 0.09500 and p = 0.09513: 0.095, 0.09504, 0.09508, 0.09511, 0.09512, 
0.095125, 0.09513, always with d = 0.5. Our estimate for the critical annihilation rate 

lb p c — U.UJOlU_ a00005 . 



3. Simulation results 



First, in figured] we present the raw simulation data of the average particle density p 
and the pair density function of time, starting from a random initial state with 

density p = 1/2 (and p p = 1/4). The particle density shows an approximately straight 
line in a double logarithmic plot, indicating power- law behavior. At the same p values, 
the pair density shows strong curvature, which indicates strong finite-time corrections 
to power-law behavior. We note that curves close to the critical point look qualitatively 
similar to the ones found in Ref. [8], although the slope is different: we find 0.33, 
whereas they reported 0.27(1). In our opinion, this discrepancy can be attributed to 
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Figure 1. The density p (black curves) and the pair density p p (gray curves) as a 
function of time in a double logarithmic plot. The lattice size is L = 2 21 , with 200 
independent runs for each value of p and 0.095 < p < 0.09513. 
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Figure 2. The effective exponent <5 c ff as obtained for both the particle density p 
(colored black) and the pair density p p (gray), as a function of t -0 - 2 . The lattice size 
is L = 2 21 , with 200 independent runs for each value of p and 0.095 < p < 0.09513. 

differences in finite-time corrections, due to details in the models. Because our data has 
very small error bars, we can use <5 e fr to investigate the "power-law" like behavior of the 
density and the pair density in more detail. 

The effective exponents 5 e g and S e s !P , describing the decay of particles and of pairs 
of particles, are plotted against t~ a in figure El with values for p close to p c . We chose 
a = 0.2, as this is a rough estimate for the exponent governing the leading finite-time 
corrections. As we will not use any extrapolation method, the choice for a does not 
affect our conclusions. 
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Figure 3. The ratio p p /p between the pair density p p and the particle density p as 
a function of t~ a 1 with a = 0.2. The precise value of a is not important here. At 
criticality the curve seems to arrive at the vertical axis at a value higher than 0, which 
means that eventually p/p p will reach a constant value. 

It is clear from the plot that 5 c s tP ^ 2<5 e ff, even at short times t, and thus we 
conclude that mean field theory does not correctly describe the TCPD model. 

The figure also shows that some curves first ascend to a value for 8 e e of mean field 
theory (1/3), then turn slightly downwards to values of around 0.32, and then ascend 
above 1/3 again. We interpret this behavior as a signature of a p- value which is close 
to, but slightly below p c , as the data roughly follow the critical curve, before reaching 
the inactive regime. Our main reason is that we find it very unlikely that the critical 
curve makes more than one bend on these time scales. Thus, we find a lower bound: 
5 < 0.32. Obviously, if our assumption is not correct this lower bound is also not correct. 
However, were this the case we believe that any attempt at numerical analysis will be 
impossible with the current simulation approach and state of software and hardware. 
Another important aspect to note is that <5 e ff and £ e ff,p are closing in on each other in a 
way very similar to that in the PCPD model [TJ. 

The ratio p p /p is plotted in figure |3] against t~ a , again with a = 0.2, for different 
values of p, including ones that are above p c . We have to be more careful here in our 
choice of a, because we are trying to make a conclusion about the extrapolated value. If 
we assume that p p /p ~ then the choice a — (3 combined with linear extrapolation 
yields the correct constant for t — > oo. On the other hand, choosing a > (3 yields an 
overestimation, whereas a < (3 gives us a prediction lower than the true value. Since 
we are not interested in the exact value of the ratio at t — > oo, but only in whether it 
is equal to or higher, we try to choose a < (3. Since (3 is unknown, we chose a value 
for a, such that the critical curves seem to arrive close to horizontally at the x-axis, 
which is indicative of a < /3. In figure [3j all curves are approaching the vertical axis 
(corresponding to t — > oo) in a way that strongly suggests that it will go to a finite 
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value. Thus, we conclude that the ratio p p /p approaches a non-zero value as t — > oo. 
4. Summary and conclusion 

We have performed extensive simulations of the one- dimensional TCPD model, using 
a highly efficient GPU-based simulation approach. We find that the TCPD model is 
not described by mean field theory, as evidenced by the convergence of the ratio of 
the pair density and the particle density p p / p to a non-zero value at criticality, instead 
of going to zero as expected by mean field theory. Given the similarities between the 
PCPD model and the TCPD model with regards to finite-time corrections, we find it 
not unlikely that they both belong to the DP universality class. Numerical evidence 
for the PCPD model belonging to the DP universality class was given by us in Ref. 
[7]. We emphasize that there is no solid numerical evidence propositioning a value of S 
for the TCPD model. However, we do find an upper bound 5 < 0.32, which excludes 
most previous literature values for this exponent in TCPD, and is also less than mean 
field. Thus a safer conclusion is that numerical data does not exclude the possibility 
that TCPD and PCPD belong to the DP universality class, and that this should still 
be considered as a serious possibility. 
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